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Abstact: For Hamiltonian systems subject to an external potential, which in the 
presence of a thermostat will reach a nonequilibrium stationary state, Dettmann 
and Morriss^^^ proved a strong conjugate pairing rule (SCPR) for pairs of Lya- 
punov exponents in the case of isokinetic (IK) stationary states which have a 
given kinetic energy. This SCPR holds for all initial phases of the system, all 
times t and all numbers of particles N. This proof was generalized by Wojtkovski 
and Liverani^"^^ to include hard interparticle potentials. A geometrical refor- 
mulation of those results is presented. The present paper proves numerically, 
Q>^ ' using periodic orbits for the Lorentz gas, that SCPR cannot hold for isoenergetic 

(IE) stationary states, which have a given total internal energy. In that case 
strong evidence is obtained for CPR to hold for large N and t, where it can be 
conjectured that the larger N, the smaller t will be. This suffices for statistical 
mechanics. 
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Q 1. Introduction. 

O ■ This paper consists of three parts: a statistical mechanical, a numerical and a mathemat- 
ical part. They are not equally accessible to readers interested in the new and developing 
^ ' interaction between dynamical systems theory and statistical mechanics of dissipative sys- 
d I tems. This reflects, in our opinion, the multidisciplinary approach involving on the one hand 
the mathematical theory of dynamical systems and on the other hand the physical theory 
of dissipative macroscopic, i.e. , statistical mechanical systems. In order to make progress 
towards a deeper understanding of the dynamical foundation of the statistical mechanics of 
dissipative macroscopic systems, such a combined approach seems unavoidable. However, it 
leads, regrettably, to difficulties in understanding the language of the different disciplines. 

Since this paper is primarily adressed to physicists, we have written the core of the paper 
in a language that should be accessible to a large portion of them. However, an interesting 
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mathematical unification and reformulation of previous results has been obtained by one of 
us (CP.). This work is of a more mathematical nature than the rest of the paper, but make 
a connection between mathematics and physics, in particular the mathematics of symplectic 
spaces and the physics of thermostatted dynamical systems. For that reason we have sum- 
marized the main results of this work in a Discussion Remark, but reserved a sketch of the 
details to an Appendix. 

The modelling of dissipative macroscopic systems by dynamical systems, i.e., by sets of 
equations of motion for the particles of the system, which mimic the behavior of such systems 
under the influence of an external field in a nonequilibrium stationary state, can be achieved 
by adding a suitable dynamical thermostatting term to the equations of motion. In fact, 
in all computer simulations the dynamical system reaches a stationary state, by adjusting a 
thermostat term in the equations of motion in such a way that the heat developed due to the 
work done on the system by the external forces is removed via effective friction terms in the 
equations of motion of the system. This approach has raised a lot of interest in the last ten 
years. In particular in 1988, Posch and Hoovert^l found a relationship between the average 
phase space contraction rate of a many particle color diffusion system under the influence 
of an external color field and the work done on the particles of the system, if the internal 
energy of the system in the stationary state was kept constant (iso-energetic (IE) thermostat). 
Another way of saying this is that in the stationary state the average phase space contraction 
rate of the system is equal to the irreversible entropy production rate in the system . This 
allowed one to express transport coeflBcients in terms of the sum of the Lyapunov exponents 
of the system and established a direct relation between dynamical systems theory in phase 
space and macroscopic physics in ordinary space. Similar relations have been found for other 
dynamical systems, including those representing viscous flow, if the stationary state is one 
with a fixed kinetic or internal energy. In the case of hard core interparticle interactions, 
when no interparticle potential energy is present, iso-energetic and iso-kinetic (IK) systems, 
which have constant kinetic energy, are the same. Most computer experiments have been 
carried out for the IK case and some for the IE case as wellt^^'t^L 

The dynamical equations describing the macroscopic systems were solved numerically, us- 
ing nonequilibrium molecular dynamics (NEMD). It was found this way that an important 
simplification occurred - for example in the computation of the transport coefficients - in 
that the sum of all Lyapunov exponents, which determines the phase space contraction rate, 
obeys a conjugate pairing rule (CPR) (see [7]). That is, the sum of every pair of conjugate 
Lyapunov exponents adds up to the same (negative) constant, which is equal to the average 
phase space contraction rate of the system in the stationary state, which in turn was equated 
to the average entropy production rate of the system. Equivalently, one can say that the 
Lyapunov spectrum of the system is symmetric. 

A mathematical proof of the CPR for the case of a possibly time-dependent, i.e., non- 
autonomous, Hamiltonian system with a constant friction added to the equations of motion, 
was given independently of the above developments in 1988 by Dressier t^l. The system she 
considered in this case was neither IK or IE. 
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For an IK Hamiltonian system consisting of N particles i = 1, with a Gaussian ther- 
mostat, which keeps the kinetic energy constant^^l, one has in three dimensions {d — 3) the 
Hamiltonian: 

N ^ 

H(P^9) = E ^ + + '^"'(^) (1-1) 
with the equations of motion 

1. Pi 

(li = — 
m 

^^ = fT\q) + fr\q))-cy{p,q)Pi 
where a{p, q) = ajxip, q) is given by 



(1.2) 



Ei Pi, [fr\q)+fr\q))] 



aiK{p,q) = (1.3) 

Pi 

where (f)^^^{q) and (p^^^{q) are the potential energy of the particles of the system and with the 
external field, respectively, where p = (pi, ...,Pn) and q = {qi, ....^q^) , pi = {px,Py,Pz) and 
qi = (q'x, <?y, qz) are three dimensional vectors, fl^^{q) — d(t)'^'^^{q)/dqi, ff^*{q) — d(t)^^^{q)/dqi, 
and (•, •) represent the usual scalar product in three dimensions, ajxip, q) is the thermostat- 
ting function (or effective friction) which will be determined in such a way that the total 
kinetic energy of the system remains constant in time. It is a mechanical representation of 
the heat exchange via the boundaries of the system representing the thermostat in a real 
system. 

Dettmann and Morriss^^^ recently proved a much stronger CPR than observed in the sta- 
tionary states of the computer simulations: their CPR holds for every initial condition 
Xq = (Pp, g^), at every time t and for any number N of particles. We refer to this Dettmann 
and Morriss pairing as a strong CPR (SCPR). Isolating two zero Lyapunov exponents, 
the remaining 6A^ — 2 local Lyapunov exponents over time t for any initial state and any 
number of particles obey the CPR. If the local Lyapunov exponents arc ordered as 
Ai > A2 > .... > A6Ar-2, then adding up the conjugate pairs Ai + A6Af-2, A2 + Agat-s, 
leads for every such pair to a constant — < a >ti where means a time time average of 
q;(p, q) over time t: — a{s)ds/t. 

It seems natural to extend the proof of Dettmann and Morriss from the IK to the IE case. 
The only difference is that a{p,q) is not given by aiK{p,q) of eq.(1.3) but by the simpler 
expression: 

E^{P^Jt^iq)) 



OiiE{p,q) = ' ^ \ ~ (1-4) 



We will call Hq the constrained part of the Hamiltonian H, i.e. 
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N ^ 

1=1 



for the IK case and 



i=l 

for the IE case. The question then arises to what extent the Dettmann and Morriss' proof of 
SCPR for the equations (1.1)-(1.3) is typical for general dissipative dynamical equations. 

The Dettmann and Morriss proof has recently been generalized by Wojtkowski and Li- 
veranit^l to include hard core contributions to the interaction potential for which 

Vi(/)*"*(Q) is not defined. 

The Dettmann and Morriss proof makes use of the observation that the CPR as stated 
above only holds in a 2dN-2 dimensional (i.e., a 2-co-dimensional) subspace of the full 2dN- 
dimensional phase space of the dynamical system. Two zero Lyapunov exponents associated 
with the conservation of kinetic energy and with the motion along the system's phase space 
trajectory have first to be eliminated, in order to allow the remaining Lyapunov exponents 
to pair to a different (non-zero) value. A similar procedure, replacing the conserved kinetic 
by the total energy does not lead to a SCPR, like that found by Dettmann and Morriss. 
In fact, computer simulations for the Lorentz gas indicate that, in general, no SCPR can 
be expected in this case. On the contrary, our computer simulations seem to indicate that 
the CPR may only hold for large systems and/or long times. The reason is not clear at the 
moment, but might be related to the special form of the friction term in eq.(1.2): —aiK{p, (l)'Pi 
which obtains both for Gaussian and Nose- Hoover thermostats!^]. For, the IK condition, 
C(ik{p, q) contains 4>^^^{q) in addition to (p^^^{q) (the latter appearing exclusively in aiE{p, q) 
(cf.eq.(1.4)). This could be responsible for the validity of the local CPR in the IK case. In fact, 
computer simulations show that the CPR continues to be obeyed to machine precision before, 
during and after a collision of one moving point particle interacting via a Weeks-Chandler- 
Andersen (WCA) potential (cf.eq.(2.1) (2.2)) with one stationary scatterer in the Lorentz 
gas IK case, but is violated considerably in the IE case. Although it is unclear at present 
why the aiK{p,q) is able to assure the CPR, the fact that it contains 0*"'*(q') = 0^^^(9), 
the same potential that causes the scattering, could well be relevant in this context. On the 
other hand, it is very hard to see how ajEip, q), which only contains the external potential, 
could physically force a CPR during a collision governed by 0*"'*(q'). 

Strictly speaking, when more than one moving particle is present in the system, the proper 
thermostatting condition involves, in general, the peculiar rather than the laboratory veloc- 
ities, that appear in the equations (1.1)- (1.4). Especially in the case of viscous flow it is 
important to identify the relevant particle velocities and the temperature of the system with 
the peculiar rather than the actual velocities of the particles, where the former refer to the 
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velocity of each particle relative to the average local flow velocity of the fluid at the position 
of the particle (cf.eqs(4.1)(4.2)). They will be used in Section 4. 
Our main results can be summarized as follows. 

1. For a Lorentz gas, where a single point particle moving in a regular triangular lattice of 
flxed spherical WCA scatterers under the influence of a (color) fleld subject to a Gaussian 
thermostat to attain a stationary state, very carefully controlled numerical simulations of 
period 3 to period 8 periodic orbits, have excluded the validity of strong SCPR for an IE 
constraint, while being consistent with an IK constraint (for all t and a;o). 

2. On the other hand both IE and the use of IK with the proper thermostatting condition 
with peculiar, rather than laboratory velocities, for 9, 16, 32 particles moving through 
WCA scatterers, strongly suggested the validity of the CPR for sufficiently large N and 
possibly for each t, with N large enough. In fact, at present it is unclear whether AT — > oo 
would be sufficient for the CPR to hold in general, i.e., for all t, or whether a large t 
condition, i.e., t — > oo, should be imposed in addition. 

The "mechanization" of the thermostats, as shown in the friction terms ~ api in the above 
equations, was developed in numerical simulations in the beginning of the 1980's'^L This 
extended to heat exchange the "mechanization" of thermal external forces - as occur in viscous 
or thermal phenomena - already introduced in linear response theory in the middle 1950'st^L 
Together they form a basis for a connection between statistical mechanical and dynamical 
systems. To what extent these "mechanized" equations of motion adequately represent the 
real physical systems and the application of dynamical systems theory considerations to them 
reflect properties of real systems, remains, at present, an open question. 

For large "mechanized" systems a number of properties of real systems, such as the positivity 
of the entropy productiont^*^! and the Onsager reciprocal relations!^ ^1 have recently been 
established. For a further discussion of this point we refer to a forthcoming papert^L 

If the generic case for the CPR to hold is large A^ and/or large t, then it is in general a 
statistical mechanical property in a nonequilibrium stationary state of a system, rather than 
a dynamical property, where it would hold for all N and/or all t. 

The outline of the paper is as follows. 

In section 2 various attempts are described to prove SCPR for the IE one particle Lorentz 
model, where one particle moves through scatterers placed on a regular triangular lattice. 
In section 3, very accurate controlled numerical results for a one particle regular IE Lorentz 
model are presented, which contain apparently incontr overt able evidence of a violation of 
the SCPR for certain periodic orbits (see Table 2), contrary to the IK Lorentz model, as 
proved by Dettmann and Morriss. In section 4, numerical evidence for an aymptotic CPR for 
large N for IE and for a different thermostat (IPK), that makes sense only for many particle 
systems and in which peculiar rather than laboratory velocities appear, is presented. Section 
5 summarizes some of the results and their implications obtained in the paper as well as a 
reformulation of the strong pairing rule condition. Appendix Al discusses the difficulties one 
encounters, if one tries to map an IE system into a Hamiltonian system, like is done in [12] 
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for IK systems. Finally Appendix A2 gives a sketch of the proof of the reformulated pairing 
rules condition in a geometrical form. 

2. Numerical results: the Lorentz model 

In this and the following sections we report the numerical simulations we have done on the 
pairing problem. We focus our attention on the simplest system: a periodic Lorentz model 
with one or more moving particles through periodically arranged scatterers. The results on 
a single moving particle are reported in this and the next section, those on many moving 
particles in sec. 4. 

The Lorentz model we consider in this section then consists of one particle moving in a 
regular triangular array of scatterers. The particle interacts with the scatterers via a Weeks- 
Chandler-Andersen (WCA) potential and is subject to an external field and an IK or IE 
thermostat. The equations of motion are thus given by eq.(1.2) where N = 1 and, writing 
direcly q, p instead of q = (qi), p = (pi), we have (j)^^^{q) — {q,E), with E a fixed three 
dimensional vector chosen to be ^ = (1.1, 0.1, 0.1), 0*''*(g) = 0^^^(g) where 0^^^(g) is the 
(finite range) WCA potential given by 

^WCA^^^ J2 ^{q-q,) (2.1) 

with 

0(^ = / MP ~ W + ll^"!'^ < ^ (2.2) 
1 otherwise 

and qn the centers of the scatterers on the regular triangular lattice. 

This system can be considered as a model for electro-conductivity in which the moving 
particle represents an electron and the scatterer represents the ion of the metal. In the 
present context it is usually used as a model for color conductivity, where the moving particle 
has a color charge rather than an electrical charge^^l, so that (electrical) interactions between 
particles, when more than one is present, are avoided. We fix the units in such a way that 
the mass of the particle is equal to 1, its charge is also equal to 1. The scatterers have a 
radius of 2^/^, a distance of 2.5 and the time the particle needs to cover a distance 1 in free 
space, i.e. when not interacting with a scatterer, is y/2. The last condition just means that 
the constant energy (kinetic or total) of the particle is fixed to 1. From a mathematical point 
of view the scatterers are needed to introduce a strong chaos in the system. In general no 
analytic results are known for this system with respect to its ergodic or chaotic properties. The 
closest system for which ergodicity and the existence and uniqueness of a limiting distribution 
(SRB) for the Liouville measure has been proved is the case in which the WCA scatterers are 
substituted by hard scatterers and the space dimension is 2 (see [13]). 

The number of degrees of freedom of the system is three. This means that we have two 
zero Lyapunov exponents and only 4 non trivial Lyapunov exponents: Ai > A2 > A3 > A4, of 
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which Ai and A2 are positive, A3 and A4 negative, if E is not too large^^^l't^'^l If we introduce 
the pairing error by: 

= Ai - A2 - A3 + A4 (2.3) 

we have for a system for which CPR holds that Pe = 0. This imphes that we can use Pe as a 
measure of the discrepancy of the Lyapunov exponents from pairing. 

Strictly speaking Lyapunov exponents are defined only in the limit when the time length of 
a trajectory goes to infinity. Moreover they are only defined almost everywhere with respect 
to some ergodic measure^^^^. In what follows we will use the term local Lyapunov exponent 
to refer to a definition of a Lyapunov exponent for a single trajectory for a finite time. The 
particular definition needed will be clear from the context. We will drop the term "local" 
when no confusion can arise. 

In our calculations for the one particle Lorentz model, we focussed our attention on three 
main points: 

A) asymptotic pairing; 

B) possibility of generalizing the DM scheme; 

C) periodic orbits (see section 3); 

A. Asymptotic pairing. 

To compute the asymptotic Lyapunov exponents, we calculated the local Lyapunov expo- 
nents for as long a time as possible, typically for about 50,000 time units. With our units the 
time unit equalled the distance unit. We use the standard algorithm'-*^^! consisting of evolving 
a basic set of unit vectors in the tangent space with the tangent flow and orthonormalizing 
them from time to time (see [14] for a theoretical justification of this algorithm). It is very 
hard to give a reliable error estimate on this kind of computation. In fact there are two main 
sources of error: 

1) the intrinsic instability of the dynamical system due to the presence of two positive Lya- 
punov exponents (quite clear from the simulation) prevents us from following a trajectory 
of the real system for a long time. In fact even using a very accurate integration method 
(we used an adaptive step size 4*'* order Runge-Kutta (RK) integrator with an error con- 
trol of lO"-'^-'^ at each time step) any truncation error will be exponentially increased. The 
typical argument to solve this problem is to consider the truncation error as a random 
perturbation and use the "shadowing lemma" to argue that there must be a trajectory of 
the real system near by the simulated one^^^^. This has, in our opinion, no real mathe- 
mathical basis. So we will proceed assuming that the Lyapunov exponents of systems like 
(2.1) behave continuously with respect to small pcrturbation-like truncation errors, even 
if we cannot prove this. We refer to [17] for further discussion. 

2) the errors introduced by the numerical round-off in the evaluation of the tangent fiow, 
in orthonormalizing the evolved basis of tangent vectors and the impossibility to do the 
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infinite time limit on a computer. To deal with this we proceeded in two different ways. 
The first one was the typical one (see e.g. [18]) and consisted in considering as error the 
variance of the computed exponent on the last part of the trajectory. The other way was 
to run more than one simulation, starting from different randomly chosen initial points. 
We got essentially the same results from both these methods. 

To have a clearer idea of the asymptotic values of the Lyapunov exponents, we monitor the 
6 local Lyapunov exponents as a function of time for a single trajectory. Using the algorithm 
of [15], we then have to follow these 6 local Lyapunov exponents. In all cases, we observe that 
after a short transient time two of them decrease to quite regularly; these two are related to 
the conservation of energy and the direction of motion. It is natural to consider the remaining 
four Lyapunov exponents as the non trivial and possibly pairing ones and use the first two 
(which we will call the zero local Lyapunov exponents) as a check on the precision reached in 
approaching the infinite time limit. In fact, only when these two zero Lyapunov exponents, 
Ao and Aqo, have become approximately equal to zero, can CPR-like behavior be expected 
to hold for the other four asymptotic- like Lyapunov exponents. We ran simulations both 
with IK and IE thermostats. Clearly we expect the pairing error for the system with the 
IK thermostat to go to zero and this can be used as a further check of the precision of the 
algorithm. The results are reported in Fig. 1. 

We see that the three lines in Fig. 1(a) for the IK system approach zero regularly. They 
can easily be described by a C/t behavior with C = 1.5 for the pairing error, C = 6.0 for the 
largest zero exponent Aq and C = for the smallest exponent Aqo- Moreover, pairing is very 
well satisfied, the pairing error being always smaller than the largest zero exponent Aq. All 
this is clearly in good agreement with [1] and reflects the fact that, if one had used a proper 
Poincare section on the energy surface and the standard definition of Lyapunov exponents 
(see [14]), one would have found an identically vanishing pairing error. In this respect, it is 
interesting to observe that the values and the fluctuations in time of the pairing error are 
always much smaller than the fluctuations of the largest zero Lyapunov exponent Aq (see Fig. 
1(c)), again indicating a strong cancellation between the exponents. 

In the IE case (Fig. lb and Id) we observe again a C/t behavior for the two zero Lyapunov 
exponents, with C = 6.0 and C = respectively. However, the pairing error does not show 
any clear behavior and strongly fluctuates, even after a very long time. Nevertheless, it must 
be recognized that its value is very small and its fluctuation is of the same order of magnitude 
as its value. We think that it is very hard to decide on the basis of these data if such a small 
value is a numerical artifact or a real effect, but we think that the fact that the fluctuations 
of the pairing error are of the same order of magnitude as those of the maximum Lyapunov 
exponent, suggest that, even if pairing is obtained, it could not be SCPR in the strong sense 
of [1], 

The results after 50,000 unit of time, with the errors computed from the variance of the last 
5,000 time units, are reported in the following table that summarizes the above discussions. 
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(a) 



(b) 



pairing error - 
first zero exponent - 
second zero exponent 



5000 1 0000 1 5000 20000 25000 30000 35000 40000 45000 50000 
t 



(c) 




painng error - 
first zero exponent - 
second zero exponent - 



5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 



(d) 



maximum exponent - 



5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 
► t 



1.08 
1.07 



maximum exponent - 



5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 



Fig. 1: Behavior, as functions of time, of (a) pairing error, Aq and Aqo as a function of time for IK Lorentz system 
(b) same for IE Lorentz system (c) maximum exponent Xmax for IK Lorentz system (d) same for IE Lorentz system 





Pe 


Ao 


Aoo 


^max 


IE 


9 • 10-^ ± 3 • 10-^ 


1 • 10-^ ± 2 • 10-'"' 


1 • 10-'* ± 9 • 10-'* 


1.0090 ± 2 • 10-^ 


IK 


4-10-^±l-10-^ 


1 • 10"^ ± 1 • 10"^ 


2- 10-6 ±3- 10-6 


1.0248 ±4- 10-^ 



Table 1: Values for pairing error (pe) zero Lyapunov exponents (Ao and Aqo) and maximum Lyapunov exponent 
i^max) for IK and IE systems after 50,000 time units 



To double check these results we also ran a simulation with 10 different initial conditions, 
both for the IK and the IE system. The results are given in Fig. 2(a),(b) where we show 
figures of the behavior of the pairing error pe for a simulation with ten initial points, plotting 
the mean value of pe over the ten points, as well as the maximum and the minimum values 
observed. We also present the behavior of the maximum Lyapunov exponent in Fig. 2(c),(d) 
to compare the error of the pairing error with that of the maximum exponent. These results 
confirm the previous discussion and the errors evaluated with this method are of the same 
order as those previously reported. Namely: for the IK systems the difference between the 
maximum and minimum values of the pairing errors is much smaller than those for the 
maximum Lyapunov exponent, while this is not the case for the IE system. It is interesting 
to note that after a long time the pairing errors for all the 10 points of the IE simulations are 
positive. Although not conclusive we do think that this can be considered as strong evidence 
of non pairing. 
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(a) 



(b) 



I 



mean value - 
minimum value - 
maximum value - 



5000 1 0000 1 5000 20000 25000 30000 35000 40000 45000 50000 
t 



(c) 



0.003 
0.0025 
0.002 
0.0015 
I 0.001 
0.0005 


-0.0005 



-0.001 



mean value 
minimum value 
maximum value 




5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 



(d) 



1.07 
1.06 
1.05 

t 1.03 
1.02 
1.01 



mean vaiue - 
minimum vaiue - 
maximum vaiue - 



1.08 
1.07 

1.06 

I 1.04 
1.03 
1.02 
1.01 



mean value - 
minimum value - 
maximum value - 




5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 



5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 



Fig. 2: Behavior of the (a) mean, maximum and minimum values of the pairing error for 10 initial conditions as a 
function of time for IK system (b) same for IE system (c) mean, maximum and minimum values of the maximum 
Lyapunov exponent for 10 initial conditions as a function of time for IK system (d) same for IE system 



B. Generalizing the DM scheme. 



In [1] Dettmann and Morriss consider the flow ipt on R^'^ generated by the Hamiltonian IK 
system (1.1)-(1.2). 

The main point in the proof of strong pairing is to split off the two zero Lyapunov exponent 
associated with the energy conservation and flow direction from the other exponents and then 
consider the pairing of the other Lyapunov exponents in an appropriate subspace of the full 
phase space. To eliminate the exponent associated with the energy it is enough to restrict the 
flow to the energy surface. To get rid of the other exponent we can use a Poincare section like 
procedure (see Appendix A2, and flgure therein, for a more detailed discussion). This means 
that given a trajectory, its initial point x(0) = (p(0),^(0)), a later point x(t) = {p{t),q{t)) 
and two corresponding planes V^(o) and V^(t) (see Fig. 10), one passing through x(0) and the 
other through x(t), both transversal (i.e. not tangent) to the trajectory, we can consider the 
map $t induced by the flow (pt that associates to a point r on V^(o) the point $t(r) on V^(t) 
which is the flrst intersection between the trajectory starting at r and the plane • If we 
properly choose a plane V^x(t) for any t, the above procedure defines a new flow $t whose 
Lyapunov spectrum is identical to the one of (ft but does not contain the zero exponents 
associated with the flow direction. This is the scheme followed in [1] where it is shown that 
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if one chooses as V^(^) the plane orthogonal to the vector Jeo(t) with J = ^ ^ J the 

usual symplectic matrix and eo{t) = {p{t),0) (we will call this choice of Vt the DM subspace), 
then the eigenvalues of the matrix Lf'Lt are paired, where the tangent flow Lt = D^t is the 
differential of (see eq.(A2.2)). From this, the pairing of the Lyapunov exponents for (pt 
follows immediately. Observe that Jeo is a vector field defined on the energy surface with no 
reference to the particular trajectory. We will call any vector field g, like Jeo above, such 
that the eigenvalues of Ll^L^, given by the above construction, are paired an a-symplectic 
field. Observe that Lf satisfies an equation of the form: 



Lt = T{x{t))Lt (2.4) 

where T(x) is called the infinitesimal generator of the tangent flow Lt. In [1] it is shown that 
the matrix Lf'Lt has paired eigenvalue if T(x) satisfies JT(x) + T{x.Y'^J = — Q!(x) J for a 
suitable a{x). If this last condition is satisfied we call Lt and the dynamics that generate Lt 
CK-symplectic and T(x) infinitesimally a-symplectic. 

An interesting property of the Lorentz model (with not too high scatterer densities, as 
considered here) is that the particle spends a significant part of its time in regions where no 
potential is present, i.e. between two successive interactions with the scatterers. In these 
regions the IE system is identical to an IK system. It is thus evident that if we consider a 
trajectory for a time small enough that the particle never interacts with a scatterer, we will 
find that the local exponents in the DM subspace are paired. Here by local exponents we 
mean those computed as the eigenvalues of ^ log L^'^L^, where Lt is the differential of the flow 
(see [14]). A natural question is now that if one considers a trajectory which starts in free 
space, goes through a collision and then arrives again in free space, can one flnd again pairing, 
using the DM subspace? This question is interesting because if there is a a-symplectic fleld 
for the IE system, this must give pairing for short (never colliding) trajectories as well as for 
long trajectories. If no pairing occurs in the DM subspace for long trajectories with collisions, 
this implies that, even if there is a cc-symplectic fleld for the IE system, there must be either 
another tranversal symplectic fleld for the IK system, (at least for the portion of phase space 
outside the scatterer potentials, i.e. in free space) or, alternatively, no DM subspace within 
which pairing exists for the IE system. 

It is not difficult to investigate this question numerically if we consider only a one collision 
trajectory and we avoid too singular collisions. For, in that case we have only to follow the 
trajectory for a short time and the diagonalization of the differential of the flow L can be 
done with very accurately so that we can achieve a precision on the local exponents of 10~^^ 
The result is that for the IE system after one collision there is no more pairing in the DM 
subspace. The interesting point is that the pairing error is of the same order of magnitude as, 
i.e. consistent with, the one observed asymptotically in the long time simulations of Fig. 1, 
discussed sub A. We were able to go through more than one collision but not through many, 
because the matrix elements of the tangent flow grow exponentially. Nevertheless for around 



11 



0.6 - 



/ 1 



-0.2 - 



IK 
IE 



0.5 



1.5 



t 



2.5 



3.5 



4.5 



Fig. 3: Behavior of tlie pairing error as function of time during a few interaction with the scatteres for the IK 
(identically 0) and IE system. The time which the particle enters and leaves the potential of a scatterer are indicated 
by vertical lines of two different dashing (the particle start outside all the potentials). 



three collisions, one sees that the pairing error in the DM space is always violated by more 
or less the same amount, most during collisions and decreasing during free flight (cf. Fig. 3). 

The question of the existence of a possibly different transversal symplectic field for the 
IE Lorentz system is nonetheless still open but this is clearly a question that is not easily 
investigated numerically. We chose to look at a similar but in principle simpler alternative: 
given a point x(0) = {p{0),q{0)) and a tangent vector w(0), at x(0) is it possible to find 
a vector w(t) tangent at x(t) such that the linearized flow from the orthogonal to w(0) to 
the orthogonal to w(t) has paired local Lyapunov exponents? This means computing the 
linearized flow for a time t starting from x(0) and then looking for a solution of the equation 
Pe{w{t)) = (where we consider the pairing error as a function of w(t)). Observe that this 
equation depends on four variables, while a unique condition is imposed so that if a solution 
exists it will generically not be unique. Note that the existence of w(0) and w(t) for any 
initial point x(0) and t does not imply that these vectors generate a vector field on the energy 
surface (like is the case in the DM proof, see Appendix A2). However, this would be enough 
to prove SCPR, if one also had that w(t) does not become orthogonal to x(t) 

The problem with this computation is that we can handle only short times t for the same rea- 
sons as given above, namely that the matrix elements of the tangent flow grow exponentially 
with the maximum Lyapunov exponent while we want to compute a function of them, i.e. the 
minimum eigenvalue, which decreases exponentially with the minimum Lyapunov exponent. 

In our simulations we followed the trajectory for a time 10 (i.e., for 5 to 7 collisions) comput- 
ing w{t) for t — nSt with St = 0.1 and n = 1, 2, 3, 4, .... The initial point was chosen outside 



12 



all the scattering potentials in such a way that the choice w(0) = Jeo(O) was a natural one. 
To obtain w(t), we started with JDHq (see (1.5), (1.6) and (A2.2)) looking for a solution 
with a steepest-descent method, till the pairing error changed sign and then using bisection. 
Observe that in this setting the algorithm applied to the IK system would be completely 
trivial. We also tried to find other solutions for example using as an initial attempt for w(t) 
at t the last computed w(t — 5t) at t — 5t. We note that this is also a possibly interesting 
procedure for the IK case because the result need not be Jeo(t), as would be the case for [1]. 

In fact, in all cases we were always able to find the vector w(t). That it was also possible 
for the IK system to find a w(t) different from Jeo(t) which satisfies the pairing rule (at 
least for a short time t) could lead to a different scheme to prove CPR for the IK system in 
a possibly more general fashion than in [1]. Unfortunately we were unable to prove this and 
our numerical simulations were too short to give real insight into the above question. 

3. Periodic orbits 

To obtain a more or less definitive result on the possibility of having pairing for all time and 
all initial conditions, we decided to check the pairing of the Lyapunov exponents on periodic 
orbits of the single particle periodic Lorentz gas, considering that periodic orbits solve the 
problem of obtaining truly asymptotic Lyapunov exponents, since the asymptotic Lyapunov 
exponents on a periodic orbit are just the logarithms of the modulus of the eigenvalues of the 
differential of the flow after a period, divided by the period. However, errors are introduced 
in finding periodic orbits numerically. Having a periodic orbit or an orbit we know is near 
enough to a periodic orbit, will permit us to obtain a very good estimate of the error involved 
in the computation of the asymptotic Lyapunov exponents. In fact, if the period of the orbit 
is not too long the exponential propagation of errors can be controlled with an accurate 
integration scheme (we again used a 4*^^ order RK with adaptive step size with error under 
10"-^^). Having to diagonalize a matrix only once, these operations do not increase the error 
significantly for sufficiently short periodic orbits, because, after a short time, the matrix 
elements of the differential of the flow are not yet very large. The real problem is to flnd 
the best approximate periodic orbits numerically and to estimate their deviations from real 
periodic orbits. 

To this end we ran a long time single particle trajectory and recorded the positions of the 
particle when it entered in the interaction (WCA) potential of a scatterer (we will register 
this as a collision). In this way we get a long list x„ of collision points and, flxing a period T, 
we look for all m such that |x^ — Xm+rl < e for some e (10~^ in our experiments). We then 
search for a periodic orbit near this orbit using the Newton-Ramson method. The algorithm 
stops when the distance between the initial and the T-th collision point is less than 10~^^. 

Having a small distance between the initial and the flnal point of the orbit does not auto- 
matically ensure that the orbit is near a periodic one. A good check of this is to see if the 
application of the Newton-Ramson method gives a quadratic convergence as should theoreti- 
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cally be the case^^^l. This means that iterating the method, if the error at the n-th step is e, 
it should be Ct^ at the (n + l)-th step, where C is a constant hnked to the first and second 
derivatives of the flow. Another good check is to see whether in the computed Lyapunov 
spectrum there are two exponents (which means that the spectrum of the differential of 
the flow must contain two eigenvalues 1 with suflicient precision. Both these conditions were 
always fulflUed in our computations: every time a periodic orbit was found, the algorithm 
converged quadratically and there were always two eigenvalues with logarithms of the order 
of 10-^2. We think that the error on the evaluated exponents can be conservatively estimated 
as less than 10~^. 

We found typically about 200 periodic orbits for periods ranging from 2 to 8. For all these 
periods it was possible to find periodic orbits whose exponents were not paired by an amount 
significativelly larger than the expected error. We report in the following table a number 
of relevant quantities we have observed for each period for the two periodic orbits with the 
largest pairing errors. 





pairing level 


pairing value 


relative pairing error 


Period 2 


-0.038903 


0.081320 


-2.090341 




0.005453 


-0.528786 


-0.010312 


Period 3 


-0.388106 


0.074560 


-0.192114 




0.061834 


-0.313116 


-0.197480 


Period 4 


-0.339775 


0.157961 


-0.464897 




0.144549 


-0.256445 


-0.563666 


Period 5 


-0.271768 


0.093026 


-0.342300 




0.084853 


-0.259349 


-0.327179 


Period 6 


-0.397877 


0.119123 


-0.299397 




0.097328 


-0.352683 


-0.275963 


Period 7 


-0.338662 


0.022185 


-0.065507 




0.020747 


-0.438584 


-0.047305 


Period 8 


-0.467454 


0.030299 


-0.064817 




0.020875 


-0.525150 


-0.039751 



Table 2: Pairing value —{a.iE)t, pairing error (eq.(2.3)) and relative pairing error (ratio of the pairing error to the 
pairing level, i.e. of columns 2 and 1) for each period of the 2 periodic orbits with the largest pairing errors. 

These results provide a numerical proof that SCPR for the one particle IE Lorentz model 
cannot be true, so that no proof like that of [1] is possible for this system. Nevertheless 
we note that this does not rule out the possibility of asymptotic pairing with respect to an 
ergodic measure. In fact, periodic orbits are a set of measure zero so that even if no periodic 
orbit has paired Lyapunov exponents, this does not imply that the asymptotic Lyapunov 
exponents for T ^ oo are not paired. 

We note that the behavior of periodic orbits can also lead to a better understanding of 
the general behavior of the system. First, there is apparently a large difference between 
the behavior of orbits with even and odd periods (even or odd refers here to the number of 
collisions on an orbit). In fact, while almost all period 3 orbits have a quite large pairing 
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error, those of even period have predominantly a zero pairing error (i.e. smaller than the 
precision of the algorithm and often actually evaluated as zero). To show this effect we 
made histograms of the probability Pr(pe) of finding an orbit with a pairing error p^, i.e. the 
number of orbits found with that pairing error divided by the total number of orbits found. 
The results are shown in Fig. 4. 
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Fig. 4: Distribution of the pairing error for orbits of period (a) 4 (b) 5 (c) 6 (d) 7. In fig. 4(a) the single point at 
about 0.65 means that about 65% of the periodic orbits with period 4 have approximately zero pairing error, while 
in Fig. 4(b) only about 1% of the periodic orbits have a pairing error near zero. The apparently layered structure 
of the histograms in Figs 4(b), (d) is due to the relatively small intervals chosen in the histograms. 



Secondly, we note that the existence of many even periodic orbits with paired Lyapunov 
exponents can perhaps be connected with some generalization of the [1] proof of pairing. We 
hope to come back to this point in a future paper. 



4. Many particle systems. 

In the previous sections we have seen that that no SCPR holds for a particular Hamiltonian 
IE system, i.e. the one particle Lorentz model. In particular, the numerical results on periodic 
orbits seem to establish that even asymptotic pairing is violated. In order to investigate the 
effect on pairing in the presence of more than one particle, in fact of many particles, we 
looked at the same model but with more than one particle moving among the scatterers. We 
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considered color diffusion in an external field, putting no interaction potential between the 
particles in order to be able to make very long and careful simulations, as well as for other 
reasons that will become clear in what follows. Observe that when a thermostat is added to 
this system of non-interacting particles, the particles are no more independent, because the 
thermostat (IK or IE) constraint introduces an indirect interaction between them. 

The equation of motions of this system are again (1.1)-(1.2) with N > 1 (we will use 
AT = 2, 3, 8, 16, 32), (jf^{q) = YliLii^i using the same notation as in sec. 2 (see comment 
above eq.(2.1)) and 0*^%) = Eii0^^^(9i) with 0^^^(^ given by eqs.(2.1)(2.2). The 
units are fixed as in sec. 2 (see comments after eq.(2.2)). 

In this setting it is possible to consider also a different kind of thermostat, which we will 
call an isopeculiar kinetic (IPK) thermostat, which is, in general, more physical, but has 
no meaning for a single particle. It consists, in constraining only the peculiar (i.e., thermal) 
kinetic energy, i.e the kinetic energy with respect to the center of mass of the moving particles. 
The equations of motion one then obtains are (i = 1, N): 

u Pi 

Qi= — 

m (4.1) 
k = /^""^(gi) + E- aiPK{p,q)~P, 

where f^^^iq) = d(f^^^{q)/partialq axid 



i=i 

is the peculiar momentum and 

aiPK{p,q) = (4.2) 

The main reason to look at a many particle system is the idea that the CPR could be 
valid only in the limit that the number of particles goes to infinity, which would suffice for 
statistical mechanics. There are two kinds of evidences for this point. The first is based on 
many numerical simulations, in which the Lyapunov spectrum of a many particle system is 
computed and a good agreement with the pairing rule is obtained, see [5] [7] [6]. However, 
considering the small deviations observed in sec. 2 and the periodic orbit results, it seems 
doubtful that those simulations can really establish pairing definitively. 

A second reason arises from an analytic argument that was first presented in [7]. We give 
here a more rigorous and extended presentation of that argument. 

Suppose we have an IE system given by a Hamiltonian of the form eq.(l.l), equations of 
motion eq.(1.2) with thermostat eq.(1.4). In Section 2, sub B, following [1], conditions for 
pairing of an IK system to hold were tied to finding a proper subspace. As was pointed 
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out before the argument for choosing the proper subspace was related to the need of taking 
out of the spectrum the two zero eigenvalues that are always present, but pair to a different 
constant (7^ 0) than the other Lyapunov exponents. However, since we have been unable to 
find the proper CPR subspace for the IE system, we were forced (as in the previous section) 
to proceed in the full phase space, replacing, in a way, the proper subspace by the condition 
that the zero Lyapunov exponents Aq ~ Aqo must be ~ in the full phase space calculations 
for CPR to hold. 

We can now consider the tangent flow associated with these equations but, instead of looking 
at it in an appropriately chosen transversal subspace, like we did in sec. 2, we look, in the 
spirit of the previous section, at the tangent flow in the full phase space. Let us call again 
(pt{x) the flow generated by eq.(1.2) and Lt = D(pt the tangent flow. Then Lt satisfles an 
equation identical to eq.(2.4), with x = ip,q): 

Lt = T(x(t))Lt (4.3) 

(we use here the same convenient notation in full phase space as we did in the DM subspace 
of sec. 2 for slightly different objects, because we think that no confusion can arise). The 
matrix T(x), the inflnitesimal generator of the tangent flow, has the simple form: 

T(x) = Th(x) -a7i,(x)/p-p0 Va7E(x) (4.4) 

where Th is the infinitesimal generator of the tangent fiow for the Hamiltonian system, i.e. 
Th = Jd^H with J the usual symplectic matrix and d^H the matrix of the second derivatives 

of H with respect to all momenta and positions. Moreover = ^ ^ || ^ , i.e. the identity 

on the p subspace and Cg> represents the tensor (or dyadic) product. It is easy to see that if 
we neglect the term Vq:/£;(x) the matrix T(x) = Tff(x) — aiE{'^)Ip satisfies the relation: 

JT(x) + T*'^(x) J = -aiE{^)J (4.5) 

where T*'"(x) is the transpose of T'(x). This follows from the fact that for a Hamiltonian 
system JTh(x) + Tff(x)*^J = and JIp + J = J. This, in turn, implies that if we 

consider the characteristic exponentst-*^^] associated with the matrix T(x), which we will call 
the mutilated spectrum of the system, they will be paired. This follows directly, for example, 
from the last part of [1], or from [7] and the pairing level will be exactly equal to — < 
aiE >t- Moreover for these exponents a SCPR holds. Note that the mutilated spectrum 
does not necessarily have two zero exponents. This creates a first difference between the 
(true) Lyapunov exponents spectrum and the (mutilated) characteristic exponents spectrum, 
but we expect that the two zero Lyapunov exponents associated with T(x) will be converted 
into a couple of paired mutilated characteristic exponents associated with T(x), one of which 
is numerically found to be zero, while the other one equals twice the pairing level — < ajE >t 
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(see Fig. 7 and discussion below it). The critical question is whether the term Vq:(x) can 
be neglected. Observe that 



daiE{^) _fr\g) 



9pi E j Pj 
9q;/e(x) 



Pi 



(4.6) 



dqi 



=0 



If we go to the thermodynamic limit, keeping the energy per particle fixed, we will have that, 
in the mean ^j'Pj — ^Wpi]], so that we can say that Vct/s is of the order of 1/N. Although 
this makes the above procedure more reasonable, it is mathematically not justified. A similar 
argument holds for IPK thermostatted systems. In fact we have in that case: 

T^^^(x) = Th(x) - aiPKi^)ip -P® Y«/px(x) (4.7) 

where Ip = Ip + 0{1/N) and Vq;7pk(x) is given by slightly more complicated expressions 
than eq.(4.6), but can still be considered as being 0{1/N). 
The l/7V-argument has three major problems: 

a) It is mathematically speaking not always true that = 0{1/N). For, in the IE case 

it can happen that 'Y^- is very small or even zero, if all the energy of the system is 
potential energy, so that the system has no kinetic energy. This implies then that the first 
term in eq.(4.6) can be very large. Of course, it will nevertheless be true that for most 
parts of phase space 'Yjpj is 0{N) and the volume of the region where it is small will 
go to zero as the number of particle goes to infinity. This effect is not present in the IPK 
system, because in the denominator of ck/pk the constrained quantity itself then appears. 

2) If we consider the two equations 

Lt = TfLt Lt = TfLt 

where f = T + 0{1/N), it is not generically true that Lf = Lt + 0{1/N) uniformly in t. 
In fact we will typically have that 



Lt = Lt + f{t)0{l/N) (4.8) 

for some f{t). The point is then to control the growth of f{t), but this is, generically, very 
hard. 

3) Even assuming that we could prove (4.8), it is easy to give an example of matrices of 
the form = + {1/N)Pn in which the spectrum of is completely difi'erent from 
that of Ln. Let L^v a Jordan block with eigenvalue a. That is {Ljs[)i^i = a, (LAr)^ = 1, 
while all other matrix elements are zero. Let {Pn)n,i = 1 while all other matrix elements 
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of Pn are zero, i.e., 



Ln = 



/a 1 
a 1 
0a 





Vo 



0\ 





a 1 
aj 



Pn = 



/O 




Vl 



0/ 



Using the properties of triangular matrices it is easy to see that the spectrum of Ln 
consists then only of the point a counted with multiplicity N, while the eigenvalues of 
Ln are given by = N~^Ui,N where Ui^N are the A^-th root of the unity, i.e. ufjy = 1. 
This implies that the spectrum of Ln converge to the circle in the complex plane of radius 
1 with center at a, when goes to infinity. Moreover, if we replace the perturbation 
0{N~^) by a perturbation 0{f{N)~^) we need essentially that limjv^oo exp{N) ~ 

the spectrum of Ln to coincide with that of Ln, i-e., any perturbation of the form 
for a fixed p will give the same result. This is clearly connected with the existence of many 
degenerate eigenvalues, but it is not clear that a hypothesis of non degeneracy would be 
enough to give convergence of the spectrum of Ln to that of Ljv. 

Since we were not able to decide analytically the above points, which are well known difficul- 
ties if one tries to prove continuity of characteristic exponents with respect to perturbations, 
we checked the behavior of the Lyapunov exponents as a function of numerically. We did 
simulations for N = 2,4, 8, 16 and 32 particles. In this case it became not feasible to run sim- 
ulations for 10 different initial conditions, because the 16 and 32 particle simulations already 
needed two weeks of CPU time. So we computed errors by only looking at the fluctuations 
in the last part of the trajectory (see sect. 2) and we always ran analogous simulations for 
the IK model to check the precision of the program. Moreover the lengths of the simulations 
were necessarily shorter than those used for the single particle system. Observe that, for the 
IE system, because of the absence of interactions among the particles, if we neglect the term 
containing the derivatives of a, the matrix T splits into a block matrix consisting on N blocks 
of dimension 6 x 6 on the principal diagonal of the matrix. This, assuming ergodicity, implies 
that the mutilated spectrum for this system consists of only 6 different characteristic expo- 
nents, each of which is N times degenerate. This allows a very easy check on the convergence 
of the spectra of the characteristic exponents associated with T, as N goes to infinity. 

Fig. 5 shows the behavior of the pairing error as a function of time, each graph reporting 
the observed pairing error for the IK, IE, IPK system at a fixed number of particles. To make 
the comparison more informative we decided to plot the behavior of the ratio He between the 
pairing error Pg and — < a >t, i.e. the value at which the exponents would pair if pairing 
were to occur. Observe that in this situation the number of pairs is no more 2, as it was in 
the one particle system. Therefore, we adopted as a definition of pairing error the variance 
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of the values observed, i.e. calling rrii — Xi + XQN-2-i we set 



(4.9) 



Observe that the behavior in time of the pairing error for the IE and IPK systems becomes 
more and more regular and similar to that of the IK system as the number of particles 
increases. In section 2 we saw that the smooth behavior of the pairing error of the form C/t 
for the IK system can be deduced from the existence of strong cancellations leading to a proof 
of SCPR. This suggests that a form of SCPR could hold when the number of particle goes to 
infinity or, at least, that the local Lyapunov exponents in a proper subspace will pair faster 
and faster in time. 
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Fig. 5: Behavior of ratio He of tlie pairing error and — < a >t as a function of time, for the IK, IE, IPK systems 
with (a) 2 particles; (b) IK, IE of (a) on a finer scale; (c) 8 particles; (d) 32 particles. 



In Fig. 6 we report the values after 5,000 time units for Hg of the three systems as a function 
of the number of particles N. 

We see that the pairing error for the IK, IE and the IPK systems is a decreasing function 
of N. Moreover the final value of the last two systems gets closer to the final value obtained 
for the IK system, where we know that the exponents should be paired. Thus it appears that 
our simulations are consistent with pairing for the IE system in the thermodynamic limit. 
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Fig. 6: Ratio lie of the pairing error Pg and — < a >t as a function of tiie number of particles for the IK, IE, IPK 
systems. The IK values provide an estimate of the computer program accuracy, since they should vanish. 



Finally, as noted before, the characteristic exponents associated with the matrix T have a 
particularly simple structure and are quite easy to evaluate. So we checked to what extent 
the Lyapunov spectrum converges to this simple spectrum. Fig. 7 reports the results for the 
IE system while Fig. 8 reports those for the IPK system. Note that in these figures only 
3N — 1 pairs of exponents are shown, because one pair can be easily deduced, i.e. it is a pair 
of zero exponents for the Lyapunov spectrum, while it is identical to the {3N — l)-th for the 
mutilated spectrum. Moreover, one would expect in general that the IE system and the IK 
system would become more and more similar as the number of particles goes to infinity. This 
is indeed seen in Fig. 9. Note that this argument does not apply to the IPK case. In fact 
this thermostat has a different and more physical nature in that the velocity of the center of 
mass is not constrained and only the energy with respect to the center of mass is fixed, so 
that we cannot expect any equivalence between the IE or IK systems and the IPK system. 

Both Figs 7 and 8 show a weak similarity of the spectrum associated with the matrix T and 
the real Lyapunov spectrum. However, the two spectra do not seem to converge to each other 
for A?" — > oo. The spectra of both models could in principle both split for A?" — > oo into six 
straight segments, an indication of which is visible in Fig. 7, at two thirds of the Lyapunov 
spectrum. Since it is rather unlikely that the spectrum associated with T will become discon- 
tinuous when A?" — > oo, it seems much more likely that the discontinuous spectrum associated 
with T will become continuous, if direct interactions between the particles are introduced. It 
is interesting to note that the two systems appear to pair to the same value and it is almost 
impossible to see from the figures whether the exponents are exactly paired or not. The 
observed pairing error of the Lyapunov spectrum is so small that any deviation of pairing 
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Fig. 7: Comparison of the Lyapunov spectrum \i of an IE system with the spectrum of characteristic exponents 
associated with T for (a) 4 particles; (b) 8 particles; (c) 16 particles; (d) 32 particles. 




could well be due to numerical errors or to too short a simulation time. The same cannot be 
said of the spectrum themselves. 

Finally we discuss Fig. 9. A generalization of the idea of the equivalence of ensembles for 
large discussed formally in a different context by a number of authors t^*^!, makes it quite 
reasonable to surmise that the IE and IK systems should become equivalent for sufficiently 
large A^. This is quite consistent with our computations. We see that the mean phase space 
contraction — < a >t is the same for the two thermostatted system, while the two spectra 
agree very well for the larger exponents, differing only slightly for the smaller exponents, in 
the latter parts of the spectra. 



5. Discussion 

1. One of the main results of this paper is the numerical proof of the non validity of SCPR for 
a one particle Lorentz gas. This is obtained through a careful evaluation of periodic orbits, 
keeping meticulously track of the computational errors. In fact, Lyapunov exponents for 
periodic orbits can be computed with very high accuracy if the periodic orbit is not too 
long. The existence of periodic orbits that do not satisfy the strong pairing rule shows 
that SCPR is not valid. We also note that periodic orbits form a set of measure zero. This 



22 



(a) 



(b) 



Lyapunov exponent 
characteristic exponent + 



4 + + 



^ * '4 



* * + + 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
► i/(2N-2) 



(c) 



Lyapunov exponent 
characteristic exponent + 



^ ^ ijs ijs ij> 



+ + + ¥ ^ 

^ >t> 'i? ^ 

+ + + + ^ 

. . « * * 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



. i/(2N-2) 



(d) 



0.5 
J^i 
I -0.5 



mm 



Lyapunov exponent 
^characteristic exponent + - 



o<> 





++ + + ++++ ^ ^ - 



0.1 0.2 0.3 0.4 



0.5 0.6 0.7 
► i/(2N-2) 



Lyapunov exponent 
characteristic exponent + - 




0.1 0.2 0.3 0.4 



0.5 0.6 0.7 
► i/(2N-2) 



Fig. 8: Comparison of the Lyapunov spectrum \i of an IPK system with the spectrum of characteristic exponents 
associated with T for (a) 4 particles (b) 8 particles (c) 16 particles (d) 32 particles. 



implies that no conclusion as to asymptotic pairing of non-periodic trajectories can be 
drawn from the above observations. Nevertheless the violation of pairing in periodic orbits 
is consistent with the results of the very long simulations for non-periodic trajectories 
of Section 2, used to compute asymptotic Lyapunov exponents. Also these simulations 
showed deviations from pairing for the IE systems, but the deviations were too small to 
identify clearly a real effect, because of the numerical errors always associated with this 
kind of simulations. 

2. We have also tried to generalize the DM scheme for proving strong pairing. This lead 
to a reformulation of the proofs in [1], [2] and [8] by one of us (CP.). With reference 
to the discussion and notation of sec. 2 sub B, we call a vector field g for the flow (pt 
defined on the energy surface a-symplectic if the eigenvalues of L*^(x)Lt(x) are paired, 
where Lf is the linearized flow for the trajectory starting at x. Then it was possible to 
give a necessary and sufficient condition on g to be a-symplectic directly in terms of the 
equations of motion . More precisely, calling the plane orthogonal to g(x) in the point 

X and T(x) = ^Lt(x) , this implies that there must exist a function q;(x), such that 

t=o 

the matrix M(x) = JT{x) + T^''{x)J satisfies: 

(w, M(x)w') = a(x)(w, w') (5.1) 



23 



(a) 



(b) 



Lyapunov IK 
Lyapunov IE 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
i/{2N-2) 



(c) 



- * * * * <^ 



* * $ 



* $ $ 



Lyapunov IK o 
Lyapunov IE + 



^ ^ ^ ^ .$> 

+ t i * * 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 
► i/(2N-2) 



(d) 



Lyapunov IK o 
Lyapunov IE + 



0.1 0.2 0.3 0.4 



0.5 0.6 0.7 0.8 0.9 
► i/(2N-2) 



1.5 
1 

0.5 
' 
-0.5 

-1 
-1.5 

-2 



1 


II 


Lyapunov IK o 
Lyapunov IE + 














1 


II 





0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
► i/(2N-2) 



Fig. 9: Comparison of the Lyapunov spectrum \i of an IE system with the characteristic spectrum of an IK system 
associated with f for (a) 4 particles; (b) 8 particles; (c) 16 particles; (d) 32 particles. 



for any vectors w, w' e V^x- For more details we refer to Appendix A2. 

3. To test the idea that for a large number of particles (i.e. in the thermodynamic limit) the 
IK and IE thermostats should give identical results we have done a number of simulations 
on a number of many particle Lorentz gases. The necessity to evaluate the full Lyapunov 
spectrum with high accuracy made it necessary to use high order integration methods 
leading to the impossibility of using a number of particles in excess of 32. In this setting 
we have considered also a different kind of thermostat (IPK) in which only the kinetic 
energy respect to the center of mass of the moving particles is constrined. The results 
of the simulations seem to be in good agreement with the hypothesis that pairing should 
hold (possibly in a strong sense) for IE and IPK systems in the thermodynamic limit. This 
hypothesis was first used heuristcally in [7]. However, as shown in section 4, it cannot 
easily be converted into a proof of CPR. Nevertheless, the numerical results obtained 
here are consistent with this hypothesis, thus suggesting an interesting property of large 
thermostatted Hamiltonian systems. 
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Al. Hamiltonian equivalence 

In [12] it was shown that an IK system for a given value of the energy can be mapped into 
a Hamiltonian system on an energy surface by a proper time rescaling and redefinition of 
the momenta. More generally it was shown in [2] that any a-symplectic dynamics can be 
mapped into a symplectic one {a — 0) one. In [12] there is no explicit derivation of the proper 
rescaling, but it is quite easy to derive this starting from a Hamiltonian system. We note 
that if such a transformation to a Hamiltonian system can be found for an IE system, this 
will quite easily lead to a proof of SCPR. Moreover we believe that, if SCPR holds for an 
IE system, some kind of mapping into a Hamiltonian system should be possible I^^^. To be 
general enough we consider a Hamiltonian H of the form eq.(l.l) with equations of motion 
given by eq. (1.2) (1.4). 

We now try to find a function t{X) of some new time variable Asuch that | '(t(A))^ + 
0*"*(g) = £ for some constant S representing the new fixed internal energy. Here the prime 
represents the derivative respect to A, i.e. g/(t(A)) = dqi{t{X)) / dX. This leads, if we restrict 
ourself to the H = E energy surface, to: 

(g) 

t W = ,. , - , (Al.l) 

^ ^ E- 0»^*(g) - 0e^*(g) ^ ^ 

Defining tt = t' {X)p{t{X)) we get 

- - {A1.2) 

, , E^^'iFr\q) 

«(7E, q) = — ^ ^ ~ 



where 



E - 0*"*(g) - 0«=^*(g) 
It is now easy to see that if we set 0*"'*(g) = we obtain 



^"'^^^ ^ = log(i? - <P^-\q)) {AlA) 



which is exactly the inverse transformation of the one in [12]. 



25 



The system (Al.3) looks very much hke an IE system. The only problem is that the force 
F^^^{q) is not a gradient unless (f)^^^{q) is a multiple of (p"^*{q) [21]. This shows that it is, 
in general, not possible to map an IE system into a Hamiltonian system, at least using the 
simple mappings discussed here. 

A2. A geometric condition for strong pairing 

In general, for the thermostatted systems we are considering, the Lyapunov spectrum con- 
tains two zero exponents, one due to the presence of a constant of the motion ( kinetic energy 
for IK, internal energy for IE and peculiar kinetic energy for IPK) and the other due to the 
motion in the time direction. To see whether the other Lyapunov exponent pair to a non 
zero value, one has to conceive a well defined procedure to split off the two zero exponents 
from the spectrum. This is done in [1] by a Poincare section like procedure. We want to 
give here a general version of this procedure, leading to an explicit algebraic condition on the 
differential of the flow for having strongly paired Lyapunov exponents. 

As a general setting we consider the phase space of our system as being a 2n — 1 sub- manifold 
S of R^'^ (the full {p, q) phase space) for some integer n defined by a particular level set of 
the constant of the motion (i.e. the set of points on which the constant of motiom takes a 
given value). The fiow itself, denoted by (pt, is generated by a vector field v, so that 

¥'t(x)=v((^t(x)) (^2.1) 

for X e S. 

We will also consider in R^^ the standard symplectic matrix — (^^j q ^ ? already defined 

in sect. 2, where / is the n x n identity matrix. It seems natural to introduce J considering 
that, both in the Hamiltonian and in the IK case pairing comes from symplectic properties 
of the flow (ft- In general, we can consider E as a 2n — 1 submanifold of some more general 
symplectic manifold Ai but we will consider only the case A4 = i?^" with the symplectic 
structure given by J. 

We now apply this to the pairing property. The Lyapunov spectrum of the flow ipt, with 
ipt defined only on E, does not contain the zero exponent due to the imposed constant of 
the motion. To eliminate also the second zero exponent in the time direction we proceed as 
follows, following [1]. At any point x e S we consider an hyperplane Fx of i?^" at x and 
uniformely transversal to Vx, i.e. the angle between and Vx is everywhere bounded away 
from 0. The intersection between and S defines a (2n — 2) -dimensional submanifold 
whose tangent space E'x is just the intersection between the tangent space of E at x, TxS, 
and Vx (here and in what follows we will identify R^^ with its tangent space at any point, 
so that TxV^ = ^)- Observe that the family can be defined as the family of hyperplanes 
orthogonal to a vector field gx on S ( in general gx need not be in TxE, it can be in i?^", 
but in what follows we will always consider gx G TxS). The family will allow us to define 
a new fiow, whose Lyapunov spectrum will not contain the second zero exponent as well. 
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Given now a point x(0) G E and its time evolved point x(t) = (/7t(x(0)) after a given time 
we can define a new map $t from a neighborhood of x(0) in nx(o) to Ox(t) as follows. 
Given a point x'(0) G fix(o) close enough to x(0), we consider the time T(x'(0),t) at wich 
the trajectory starting from x'(0) intersects ^x{t)- Here and in the following we suppress 
the dependence on x(0). Setting then ^>i(x'(0)) = </3T-(x'(o),t)(x'(0)) it is clear that ^t(x) 
can be considered as a new flow and that its Lyapunov spectrum, i.e. the eigenvalues of 
A = limt_»oo -^^og D^f'D^ti are identical to those of (pt at x(0), except that they do not 
contain the two zero exponents. Here D^t : V^(o) — ^ ^(t) is the differential of $t and is given 
in local coordinate by: 




{A2.2) 



Fig. 10: Schematic representation of the construction of the flow 4>t. Here x(t) = <pt(x(0)), x'(t) = (pt(x'(0)) and 
x'(t) = ^>t(x'(0)) e ^^(t)- The other symbols are explained in the text. 

Oberve that, calling Lt = D^t^ we have Lt = T{yi{t))Lt where T(x) is a matrix uniquely 
defined at every point x of S, independently of x(0), contrary to $t and L^. In fact we have: 

nx).. = |(..,(x)),, = . ^) (-3) 

so that, for w G one has, after straightforward algebraic manipulations: 
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T(x)w = 



V, + 



(Dr) 



w 



d d . 

^v. + -(i^c,,)w 



where Dt is the differential of r with respect to x. 
Noting that {Dipt) — (-^v) [Dipt] we obtain: 



T(x)w = 



dt 



-w 



It is now easy to show that = b + b with 



{A2A) 



(A2.5) 



b = 



g + (Dv) g 
(v, g) 



(yl2.6) 



and b||g. Moreover the component of the vector b along is unique. Thus eq.(A2.5) can 
be written as: 



Tw = {Dv)w + (b, w)v (^2.7) 

where w G and b is given by eq.(A2.6). 
It is known, see [1] for a simple proof, that if for any w, w' G i?xo 

{Ltw, JLt^') = ;Ui(w, Jw') (A2.8) 

where Ht = exp(jQ Q;(x(r))dr) for some function a{x.) defined on E, then if A is in the 
spectrum of At = Lf^Lt, so is /U^/A. This immediately implies that if limt^oo ^ log At exists 
then its eigenvalues are strongly paired in the sense of sec. 1. Moreover, for any ergodic 
measure for (ft the Lyapunov exponents of (ft are paired. Observe that in order give meaning 
to eq.(A2.8) we need that is a symplectic plane, i.e. JV^ — Vx- To achieve this we must 
set g = Jn where n is the normal vector to S at x. We will say that Lf is cu-symplectic if 
eq.(A2.8) holds. Moreover we will say that the flow (pt is a-symplectic if for any xq and for 
any t the matrix Lf is a-symplectic. We can now state the result. 

Theorem A2.1.- The flow (pt is a-symplectic if and only if for any x G S and for all 
w, w' G E-^ the skew symmetric matrix M = JT + T^^J satisfies: 

{w, Mw') = a{w, Jw') {A2.9) 
where a : T, ^ R is a suitable function and T is given by eq.(A2.7). 

Eq. (A2.9) follows immediately by differentiating eq.(A2.8) with respect to t. 
Theorem A2.1 provides a well defined algebraic condition that can be easily checked in cases 
of interest. In fact, straightforward computation shows that it is indeed satisfied in the case 



28 



of the IK dynamics studied by Dettmann and Morris in [1] as well as in the case of constant 
a, studied by Dressier in [8]. 
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